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Abstract 

We investigate the effects of relativity on the gravitational instability of finite isother- 
mal gaseous spheres. In the first part of the paper, we treat the gravitational field within 
the framework of Newtonian mechanics but we use a relativistic equation of state in 
the condition of hydrostatic equilibrium. In the second part of the paper, we study the 
full general relativistic problem for a gas described by an equation of state p = qe such 
that the pressure is proportional to the energy density ("isothermal" distribution). For 
q = 1/3, this equation of state describes the core of neutron stars. The mass-density di- 
agram displays some damped oscillations and there exists a critical value of mass-energy 
above which no equilibrium state is possible. We show analytically that the mass peaks 
are associated with new modes of instability. These results are strikingly similar to those 
obtained by Antonov [Vest. Leningr. Gos. Univ. 7, 135 (1962)] and Lynden-Bell & Wood 
(1968) for a classical isothermal gas. Our study completes the analogy between isothermal 
spheres and neutron stars investigated by Yabushita [MNRAS 167, 95 (1974)]. 

1 Introduction 

Isothermal spheres play an important role in astrophysics. They were initially introduced in 
the context of stellar structure (Chandrasekhar 1942) when composite configurations of stars 
consisting of an isothermal core and a polytropic envelope were constructed and studied. On 
larger scales, they were applied to stellar systems such as globular clusters and elliptical galaxies 
(Binney & Tremaine 1987). The age of globular clusters is such that their isothermal structure 
is due to a succession of encounters between stars that lead to an equipartition of energy, 
like in an ordinary gas. This statistical mechanics prediction has been confirmed by direct 
observations of globular clusters and is well reproduced by King's models that incorporate 
a truncation in the distribution function so as to account for tidal effects. In the case of 
elliptical galaxies (and possibly other collisionless self-gravitating systems like Dark Matter), 
the isothermal distribution is a result of a "violent relaxation" by phase mixing (Lynden-Bell 
1967). In that case, there is no segregation by mass contrary to the collisional relaxation. 
In general, violent relaxation is incomplete and the distribution function must be modified at 
high energies. However, the distribution function that prevails in the inner regions of elliptical 
galaxies is isothermal and this is an important ingredient to understand de Vaucouleurs' R 1 ^ 
law (Hjorth & Madsen 1993). Isothermal distributions are also appropriate to describe the cold 
interstellar medium where the temperature is imposed by the cosmic background radiation at 
T ~ 3K in the outer parts of galaxies, devoid of any star and heating sources (Pfenniger & 



Combes 1994, de Vega et al. 1996). They can also be of interest in cosmology to understand 
the fractal structure of the universe (Saslaw & Hamilton 1984, de Vega et al. 1998). 

The stability of isothermal gaseous spheres in Newtonian gravity has been investigated by 
Antonov (1962), Lynden-Bell & Wood (1968) and Katz (1978) by using thermodynamical argu- 
ments and topological properties of the equilibrium phase diagram by analogy with Poincare's 
theory of linear series of equilibrium. They showed in particular that, if the energy or the 
temperature are below a certain threshold, no hydrostatic equilibrium can exist. In that case, 
the system is expected to undergo a phase transition and collapse. This is called "gravothermal 
catastrophe" when the energy is kept fixed (microcanonical ensemble) and "isothermal collapse" 
when the system evolves at a fixed temperature (canonical ensemble). The thermodynamics 
of self-gravitating systems was reconsidered by Padmanabhan (1989) who calculated explicitly 
the second order variations of entropy and reduced the problem of stability to an eigenvalue 
equation. His method was applied in the canonical ensemble by Chavanis (2001) who showed in 
addition the equivalence between thermodynamical stability and dynamical stability based on 
the Navier-Stokes equations (Jeans problem). Similar studies have been performed by Semelin 
et al. (1999,2000) and de Vega & Sanchez (2001) by using field theoretical methods. As shown 
by Padmanabhan (1989) and Chavanis (2001), the problem of stability can be studied without 
approximation almost analytically (or with graphical constructions) by making use explicitly 
of the Milne variables introduced long ago in the context of stellar structure (Chandrasekhar 
1942). 

We show in the present paper that these methods naturally extend to the context of gen- 
eral relativity. Isothermal gaseous spheres (in the sense given below) have been only poorly 
studied in general relativity despite their similarity with classical isothermal spheres. The most 
extended study that we have found is the contribution of Chandrasekhar at the conference 
given in the honour of J.L. Synge in 1972. His work was completed by Yabushita (1973,1974) 
who considered the stability of a relativistic isothermal gas surrounded by an external envelope 
imposing a constant pressure (an extension of the classical Bonnor (1956) problem). In order to 
make the link with the works of Antonov (1962) and Lynden-Bell & Wood (1968) in Newtonian 
gravity, we shall consider the situation in which the volume is fixed instead of the pressure. 
This study can have direct applications to the stability of neutron stars since the equation of 
state that prevails in the central region of these highly relativistic objects is very close to the 
"isothermal" one. 

In section ^j, we consider the stability of isothermal gas spheres described by a relativistic 
equation of state but still using the framework of Newtonian gravity. This is a first attempt 
to introduce relativistic effects in the classical Antonov problem. In fact, as is well-known, 
special relativity does not modify the equation of state for a perfect gas (Chandrasekhar 1942). 
Only does it alter the onset of instability. We find that the critical energy below which no 
hydrostatic equilibrium is possible depends on a relativistic parameter /z defined as the ratio 
between the size of the domain R and a "classical" Schwarzschild radius R c J; ass = 2GM/c 2 , 
where M is the total mass of the system. In the classical limit R ^> Rg ass , we recover the 
Antonov result E c = — 0.335GM 2 / R but when relativistic corrections are included we find that 
the critical energy is increased, i.e., instability occurs sooner than in the classical case. The 
density perturbation profile that triggers the instability is calculated explicitly. For fi — > +oo, 
it presents a "core-halo" structure but relativistic effects tend to reduce the extent of the halo. 
When the system is maintained at a fixed temperature instead of a fixed energy (canonical 
description), the results are unchanged with respect to the Newtonian case. 

Of course, when the Schwarzschild radius becomes comparable to the size of the system, 
general relativistic effects must be taken into account. This problem is treated in detail in 
section [| We consider a simple equation of state p = qe (where q is a constant) which generalizes 
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the equation of state for isothermal spheres in the Newtonian context (the classical limit is 
recovered for q — > 0). Since the equations governing equilibrium have the same structure 
and the same properties as in the classical case, we shall say that the system is "isothermal" 
(following the terminology of Chandrasekhar 1972), although this equation of state does not 
correspond to thermal equilibrium in a strict sense (see section |3~2| ). True statistical equilibria 
in general relativity have been investigated by Katz et al. (1975) and they are characterized by 
a non uniform temperature (because of the gravitational redshift). However, systems described 
by the equation of state p = qe are numerous in nature and they include for example the 
important case of neutron cores which are usually modeled as an assembly of cold degenerate 
fermions for which q = 1/3. This simple isothermal equation of state is the high-density limit 
of more general equations of state usually considered for neutron stars (Oppenheimer & Volkoff 
1939, Meltzer & Thorne 1966). Quite remarkably, our simple "box" model is able to reproduce 
the main properties of these objects. This suggests that the structure of neutron stars is due 
intrinsically to their isothermal cores and not to the details of their envelopes (which are not 
well-known) . In section |3.5| , we show that the mass-density diagram displays an infinity of mass 
peaks and that isothermal spheres exist only below a limiting mass corresponding to the first 
peak. In section |3.6| , we show that the mass-radius diagram has a spiral behavior similar to the 
one observed for neutron stars. Using an equation of pulsation derived by Yabushita (1973), 
we demonstrate analytically that the series of equilibrium becomes unstable precisely at the 
point of maximum mass (like in the study of Misner & Zapolsky 1964) and that new modes of 
instability correspond to secondary mass peaks. We obtain the same stability criterion from 
energy considerations based on the binding energy E = M — Nmc 2 . Said differently, these 
results indicate that, for a fixed mass, the system becomes unstable when its radius is smaller 
than a multiple of the Schwarzschild radius, a property consistent with Chandrasekhar 's (1964) 
general theory. The perturbation profiles of density and velocity are calculated explicitly and 
expressed in terms of the Milne variables. They do not present a "core-halo" structure. All 
these properties are strikingly similar to those obtained for classical isothermal spheres in the 
canonical ensemble (see, e.g., Chavanis 2001). This completes the analogy between isothermal 
spheres and neutron stars investigated by Yabushita (1974). 



2 Antonov instability for a relativistic gas sphere 
2.1 The relativistic gas 

We consider a system of N particles, each of mass m, in gravitational interaction. We allow 
the speed of the particles to be close to the velocity of light so that special relativity must be 
taken into account. However, in this first approach, we shall treat the gravitational field within 
the framework of Newtonian mechanics. This procedure is permissible if the typical size of 
the system is much larger than the Schwarzschild radius (see, e.g, Chandrasekhar 1942). Let 
f(r,p,t) denote the distribution function of the system, i.e. /(r, p, t)d 3 rd 3 p gives the average 
mass of particles whose positions and momenta are in the cell (r, p; r + d 3 r, p + d 3 p) at time t. 
The integral of / over the momenta determines the spatial density 

(1) P = J fd 3 p } 
and the total mass is expressed as 

(2) M = [ pd 3 r. 
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In a meanfield approximation, the total energy is given by 

(3) E = [ — e d 3 rd 3 p + - f p$d 3 r = K + W, 

J m 2 J 

where K and W are the kinetic and potential energy respectively. According to the theory of 
special relativity, the energy of a particule reads 

This expression does not include the rest mass, so it reduces to the usual kinetic energy ^ in 
the Newtonian limit. All the effects of gravity are incorporated in the potential energy which 
contains the gravitational potential $ related to the density by the Newton-Poisson equation 

(5) A$ = AnGp. 

We now ask which configuration maximizes the Boltzmann entropy 

(6) S=-k [ — In —d 3 rd 3 p, 

J m m 

subject to the conservation of mass M and energy E. To that purpose, we proceed in two steps. 
In this section, we maximize S[f] at fixed E, M and p(r). This provides an optimal distribution 
function expressed in terms of p. Then, in section [2.2| , we maximize S max [p] at fixed M and E. 



This optimization problem is not trivial and will be further discussed in sections |2.4| - |2T6 . 

Since the gravitational potential can be deduced from the density by solving the Poisson 
equation fl5|), maximizing S at fixed E, M and p(r) is equivalent to maximizing S at fixed p(r) 
and K. Writing the variational principle in the form 

(7) SS - kf35K - J \(r)5pd 3 r = 0, 

we obtain the optimal distribution function 

(8) /(r,p)=A(r)e- /3e , 

which is a global entropy maximum with the previous constraints. Eq. (|j) is the relativistic 
Maxwell-Boltzmann distribution with an inverse temperature 

The Lagrange multipliers A(r) and (3 must be related to the constraints p(r) and K. As 
discussed in detail by Chandrasekhar (1942), these relations can be expressed in terms of the 
modified Bessel functions 

f+OO 

(10) K n {z)= e- zcoshe cosh{n6)d6. 

Jo 

Using Eqs. (|T]) (§) (^) and introducing the Juttner transformation -2- = sinh 6 in the integral, 
the density can be written 

4-7TT77 

(11) p(r) = —^A(r)K 2 (x)e x , 

x 
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where we have introduced the dimensionless parameter 



(12) x = (3mc 2 = —. 



This parameter quantifies the importance of relativistic effects. The classical limit corresponds 
to x — > +00 (kT <C mc 2 ) and the ultra-relativistic limit to x — > (&T ^> mc 2 ). 
The distribution function (§) can now be expressed in terms of the density as 

(13 > ^• p ) = S^7S) p(r)e "'"' 

In the classical limit (x — > +00) we recover the standard formula 

and in the ultra-relativistic limit (x — > 0), we get 

c 3 

(15) /(r ' p) = 8^T^ P(r)e ""- 

Similarly, after some elementary transformations, the kinetic energy can be expressed in 
terms of the normalized inverse temperature x by 

(16) K = T{x)Mc 2 1 with f(x) = 3K3{x i + , Kl{x) - I. 
Using the recursion formula 

In 

(17) K n ^(z) - K n+1 (z) = K n (z), 

the function J-(x) can be written in the equivalent form 

It has the asymptotic behaviors 

(19) ^(x) ~ A (x - +00), 

2x 



(20) :F(x) ~ - (x -> 0). 

x 

When substituted in Eq. (0), we recover the usual expressions for the kinetic energy K = 
^NkT in the classical limit and K = 3NkT in the ultra-relativistic limit. 

We can now express the entropy in terms of the density p and the inverse temperature x. 
Substituting the optimal distribution function (|T3|) in Eq. (Q), we get, up to an additional 
constant 



(21) S = kNQ(x) -k [ — In — d 3 r 



m m 
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where 

(22) Q(x) = x[J r (x) + l]+\nK 2 {x)-hix. 

The function Q(x) has the asymptotic behaviors 

3 

(23) Q(x) ~ — -hire (x — > +oo), 

(24) 31nx (z->0). 

The thermal contribution to the entropy in the classical limit is = -fciVlnT and in the 
ultra-relativistic limit Sth = 3kN In T. 

There exists a general relation between the derivatives of T and Q that we shall need in the 
following. Differentiating Eq. with respect to x and using the identity 



TL 

(25) = -K n ^(x) - -K n (x), 

x 

for n = 2, we find that 

(26) = xF'{x). 

2.2 First and second order variations of entropy 

In the preceding section, we have expressed the entropy and the kinetic energy in terms of the 
density p(r) and the temperature T (through the variable x). We now wish to maximize the 
entropy S[p] at fixed E and M. For convenience, we shall introduce a new variable y = T{x). 
In terms of this variable, the total energy and the entropy can be written 

(27) E = yMc 2 + p$d 3 r, 

(28) S = kNg(x{y)) - k [ — In ^d 3 r. 

J m m 

We can now determine the variations of S around a given density profile p(r). To second order 
in the expansion, we get 



(29) , s = + fs P (i + u£.)fy->L fm 

ay ay 2 2 m J \ m J m J 2p 

Using the identity (PE|), we find that 

dQ _dQ dx _ Q'(x) _ 
dy dx dy J- '(x) 

Differentiating one more time with respect to y, we obtain 

d 2 Q dx 1 



(31) 



dy 2 dy ^(x) 
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Substituting the above results in Eq. fl29|), we get 

(32) bS = kNx5y + kN-^ { ^ - - I ( 1 + In £-) bpd 3 v - - I ^d 3 r. 

J- (x) 2 m J \ m J m J 2p 

We now need to express the variation by in terms of bp. From the conservation of energy, we 
have the exact identity 



(33) 



= bE = Mc 2 by + J <Mpd 3 r + i / bpb§d 3 r. 



Substituting the foregoing expression for by from Eq. ( p3|) in Eq. (|32|) , we obtain 

(34) ~— f( 1 + ln— ) SpA -- / 

m J \ m J m J 2p 

Introducing a Lagrange multiplier a to satisfy the conservation of mass, the condition that 
S is an extremum is written (to first order) 



(35) = bS - abM 



- + — 1 + ln-^-J+a 
1 m\ m 



bpd 3 v. 



This condition must be satisfied for any variations bp. This yields the Boltzmann distribution 

(36) p = Ae-^, 

like for a classical gas. The condition that the critical point (^) is an entropy maximum 
requires that 

(37) b 2 S=-— [ bpb$d 3 r + ~ l — ( f$5pd 3 r) 2 -- [ ^d 3 r < 
v ; 2T J H 2Mmc 4 r(x)\J H ) m J 2p 

for any variation bp that conserves mass to first order. 

2.3 The Virial theorem for a relativistic gas 

In this section, we derive the form of the Virial theorem appropriate to an isothermal gas 
described by a relativistic equation of state. Quite generally, the potential energy of a self- 
gravitating system can be expressed in the form (Binney & Tremaine 1987) 



(38) W = - j prV$(fr. 

If the system is in hydrostatic equilibrium, then 

(39) = -pV$. 

Substituting this identity in Eq. (^) and integrating by parts, we get 

(40) W = rVpdh = i prdS - 3 / pd 3 r. 



If the pressure p b on the boundary of the system is uniform, we can write 
(41) (p prdS = Pb f rdS = p b V.rd 3 r = 3p b V, 



where V is the total volume of the system. Therefore, for any system in hydrostatic equilibrium, 
one has 

(42) W = 3p b V - 3 / pd 3 r, 



which can be considered as the general form of the Virial theorem for self-gravitating systems. 
Now, the pressure of an ideal gas can be expressed as (Chandrasekhar 1942) 

1 f /_> - 



(43) p=- -P^d'p. 

3 J m op 

For a relativistic gas described by the distribution function we get 

4tt . _ r +o ° d 



(44) p =-^wy„ 5^> 3 *- 

Integrating by parts, we obtain 

(45) p = ±- = ^kT. 

mp m 

Therefore, the equation of state for a (non quantum) relativistic gas is the same as for a classical 
gas. Written in the form 

(46) W + 3NkT = 3p b V, 



the Virial theorem ([42]) also has the same form as for an isothermal classical gas. However, the 
Virial theorem is usually expressed in terms of the kinetic energy instead of the temperature. 
Therefore, the appropriate form of the relativistic Virial theorem reads 

(47) W + = 3 Pi V. 

X \ Mc 2 I 

In the classical limit [x — * +oo), it reduces to the well-known formula 

(48) W + 2K = 3p b V, 
and in the ultra- relativistic limit (x — > 0), we get 

(49) W + K = 3p b V. 

It should be emphasized that Eq. ( f4"T| ) is only valid for a relativistic gas in thermal equilibrium. 
2.4 The equilibrium phase diagram 

Since the equation of state for a relativistic gas is the same as for a classical gas, the equilibrium 
configurations of such systems correspond to the isothermal gas spheres described extensively 
in the monograph of Chandrasekhar (1942). Only the onset of the gravitational instability will 
be modified by relativistic effects. 
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For non rotating systems, the equilibrium states are expected to be spherically symmetric. 
In that case, the Poisson equation (|5|) together with the Boltzmann distribution ( |3"6| ) yield the 
second order differential equation 

(50) -A{r 2d f)=^GAe-^. 

r 2 dr \ dr J 

This equation can also be deduced from the condition of hydrostatic equilibrium ( p9|) when 
the pressure is related to the density according to the equation of state ([45]). It is well-known 
that the density profile of such isothermal configuations behaves like p ~ r~ 2 at large distances 
so their total mass is infinite. Following Antonov (1962), we shall avoid this infinite mass 
problem by confining artificially the system within a spherical box of radius R. It is only 
under this simplifying assumption that a rigorous thermodynamics of self-gravitating systems 
can be carried out (see, e.g., Padmanabhan 1990, Chavanis et al. 2001). This procedure is 
justified physically by the realization that a distribution of matter never extends to infinity so 
R represents an upper cut-off at which other processes intervene to limitate the spatial extent 
of the system. Of course, different cut-offs are possible but fixing the volume is consistent with 
the traditional viewpoint of statistical mechanics and it is sufficient to capture the essential 
physics of the system (see the different comparisons of truncated models performed by Katz 
(1980)). 

We now wish to determine the equilibrium phase diagram of a relativistic isothermal gas. 
To that purpose, we introduce the function ip = /3m($ — $ ) where $o is the gravitational 



potential at r = 0. Then, the density field ([36|) can be written 
(51) p = p e^, 

where p is the central density. Introducing the notation £ = (AnGPmpoY^r, the Boltzmann- 
Poisson equation (|50| ) reduces to the standard Emden form 



< 52 > 

Eq. d52|) has a simple analytic solution, the singular sphere 



(53) e"* = ^, 

whose central density is infinite. The regular solutions of Eq. fl52"|) must satisfy the boundary 
conditions 

(54) ^(0) = ^(0) = 0, 

at the center of the sphere. These regular solutions must be computed numerically. In the case 
of bounded isothermal systems, we must stop the integration at the normalized box radius 

(55) a = (47rG(3mp ) 1/2 R. 

We shall now express the parameter a in terms of the temperature and the energy. According 
to the Poisson equation @, we have 

(56) GM = f R AnGpr 2 dr = f ^- (r 2 ^-) dr = (r 2 ^- 

Jo Jo dr \ dr J \ dr j r=R 
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which is just a particular case of the Gauss theorem. Introducing the dimensionless variables 
defined previously, we get 

(K>7\ _(3GMm 

(57) r] = — - — = ai[> (a). 

The relation between a and the normalized temperature 77 is not affected by special relativity. 

For the total energy, using the Virial theorem ([46]) and the expression fli~6|) for the kinetic 
term, we have 

37V 

(58) E = K + W = T(x)Mc 2 - — + 3p(R)V. 
Now, the pressure at the boundary of the domain can be written 

(59) m = m = e^. 

mp mp 

Expressing the central density in terms of a, using Eq. (|55|), we have equivalently 

2 

The total energy therefore reads 

ER Rc 2 3 
( } = ~~GM 2 ~ ~~GM + ^/{a) ~ 

where we have used Eq. (Elf) to eliminate the temperature in the last two terms. 
It will be convenient in the following to introduce the parameter 

which is twice the ratio between the system radius R and the "classical" Schwarzschild radius 

2GM 



(63) R 



class _ 

S ~~ 



Clearly, our semi-relativistic treatment (ignoring general relativity) is only valid for /x ^> 1. 
However, in our rather formal analysis, we shall treat \x as a free parameter varying in the 
range < /x < +00. The relativistic parameters /x and x are related by 

(64) x = /IT]. 

Therefore, in accordance with Eq. (|61~D and (|64D , the relation between the parameter a and 
the normalized energy A takes the form 

ER 3 e -^(«) 

(65 > A S -GAP = -<^ a <^ +^ m ~ WW - 

In the classical limit (/x — > +00), we recover the result of Lynden-Bell & Wood (1968) 

3 e'^ 



(66) A 
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Figure V. Equilibrium phase diagram for isothermal gas spheres described by a relativistic 
equation of state. Relativistic effects shift the onset of instability to larger energies. 

and in the formal limit fi — > 0, we get 

e -<M Q ) 

< 67 > A = -^- 

On Fig. [1|, we have represented the equilibrium phase diagram (3 — E for different values 
of the relativistic parameter fi. We see that the effect of relativity is to shift the spiral to the 
left. Therefore, the "gravothermal catastrophe" corresponding to the absence of equilibrium 
below a critical energy occurs sooner than in the Newtonian case. The critical energy is ploted 
as a function of the relativistic parameter fi on Fig. ^. For fi — > +oo, we recover the classical 
result of Antonov A c (+oo) = 0.335 and, in the formal limit fi — > 0, we get A c (0) = —0.345. 
Clearly, the spiral is not destroyed by relativistic effects. Only is its shape slightly modified: 
the relativistic spirals are more "stretched" than the classical one. Note that the maximum 
value of f] is independant on fi and is equal to its classical value r] c = 2.52. On the other hand, 
substituting the expression ( |5~3"D for the singular sphere in Eqs. (|57D-(|63D, we find that the 
center of the spiral is defined by the equations 

(68) Vs = 2, 

(69) A,(/i) = l-^(2fi). 

In the limit fi — > +oo, A s (+oo) = - and in the limit fx — > 0, A s (0) = — \. 

The spiral is parametrized by the normalized box radius a that goes from (ordinary gas) 
to +oo (singular sphere) when we spiral inwards. If one prefers, we can use a parametrization 
in terms of the density contrast 

(70) K = Po = e^ a) 
1 } p{R) 

that goes from 1 to +oo. On Fig. |], we plot the critical density contrast (corresponding to 
A c ) as a function of the relativistic parameter fi. For fi +oo, we recover the classical value 
ft c (+oo) = 709 (and a c (+oo) = 34.4). For fi = 0, we get K c (0) = 132 (and a c (0) = 16.0). It 
is found that instability occurs for smaller density contrasts when relativity is accounted for. 
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Figure 2: Critical energy A c as a function of the relativistic parameter /i. 



600 



3 400 





1 1 1 1 


' I ' I 

/ R o =709 






/ a c =34.4 






/ (classical) 




R c =132 / 






a c =16.0 












I , I 





-4 -2 2 4 

logfoi) 

ure 3: Critical density contrast 1Z C as a function of the relativistic parameter 
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2.5 The Milne variables 



It will be convenient in the following to introduce the Milne variables (u, v) defined by (Chan- 
drasekhar 1942) 

(71) u = — ■ — , and v = Eib' . 

ip' 

Taking the logarithmic derivative of u and v with respect to £ and using Eq. (|52|) , we get 

(72) ^H< 3 — )• 

Taking the ratio of these equations, we find that the variables u and t> are related to each other 
by a first order differential equation 

. u dv u — 1 

(74) — — = . 

v du u + v — 3 

The solution curve in the (u, v) plane is well-known and is represented on Fig. |j. Its striking 
oscillating behavior has been described by a number of authors (see in particular Chandrasekhar 
1942). We refer to Padmanbhan (1989) and Chavanis (2001) for the description of its main 
characteristics in connexion with the present work. 

It turns out that the normalized temperature and the normalized energy can be expressed 
very simply in terms of the values of u and v at the normalized box radius a. Indeed, writing 
u = u(a) and v = v(a) and using Eqs. fl57j ) (§51), we get 

(75) rj = v , 

(76) A = -fiFQiVo) + - - ^. 

^0 ^0 

The intersections between the curves defined by Eq. ([75])- (|7fi|) and the spiral in the (it, v) 
plane determine the values of a corresponding to a given temperature or energy. Considering 
Eq. (fT5|), we find that there is no intersection for 77 = ^ ^ > v max = 2.52. In the canonical 
ensemble, a gaseous sphere is expected to collapse below a critical temperature kT c = ^^r - 
This classical result is not altered by special relativity. Considering now the microcanonical 
ensemble, we first note that Eq. ( [75| ) can be rewritten 

(77) m = 3 - nvofi/jvo) - Av . 

In the classical limit (ji — > +00) it reduces to the straight line found by Padmanbhan (1989) 

(78) u = - - Av , 

and in the limit fi — > 0, we find another straight line 

(79) u = -Av . 

The curve ( [7^) is ploted on Fig. [| for a fixed value of fi and for different values of A. For 
A > A c (/i) there is no intersection, for A = A c (/x) the curve ([77|) is tangent to the spiral and 
for A = A c (/x) there are one or several intersections. We recover therefore by this graphical 
construction the existence of a critical energy below which no hydrostatic equilibrium can exist 
for isothermal spheres. On Fig. [|, we plot the same diagram as Fig. |] but for different values 
of n and, in each case, for the critical energy A c (ji). It confirms that instability occurs sooner 
when relativistic effects are taken into account. 
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Figure 4: The (u, v) plane. All isothermal spheres must necessarily lie on the spiral. There 
exists solutions in the canonical ensemble only for rj < 2.52. In the microcanonical ensemble, 
the critical energy depends on the relativistic parameter \i. For \i — 1, there exists solutions 
only for A < A c = 0.06 
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Figure 5: Same as Fig. |]but for different values of [i and, in each case, for the critical parameter 
A c (/i) above which there is no equilibrium solution. 
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2.6 The condition of t her mo dynamical stability 



We now address the question of thermodynamical stability. Let us recall that the isothermal 
spheres lying on the spiral are critical points of entropy but they are not necessarily entropy 
maxima. To determine whether they are local entropy maxima or saddle points, we must 
examine the sign of the second variations of entropy given by Eq. (|5T|) . In fact, this condition 
of stability has the same form as in the classical case provided that we make the substitution 



(80) 



1 



k 



1 



3NkT 2 



2Mmc 4 F(x)' 



Therefore, the analysis of Padmanbhan (1989) for the classical Antonov instability can be ex- 
tended straightforwardly. Introducing the mass perturbation q(r) = 6M(r) = J r 47rr 2 5p{r')dr' 
within the sphere of radius r such that 



(81) 



Sp 



1 dq 



Attt 2 dr ' 

the second variations of entropy can be put in a quadratic form 



(82) 
with 

(83) K(r,r') 



k 



5 2 S 



1 



R rR 



drdr'q(r)K(r, r')g(r') 



2Mmc 4 Fix) dr dr ' + 2 ^ 



" G 


k d 


f 1 


dY 


Tr 2 


m dr 


\Arrpr 2 


dr J 



Clearly, the conservation of mass imposes the boundary conditions q(0) = q(R) = 0. The 
problem of stability can therefore be reduced to the study of the eigenvalue equation 



(84) 



R 



K(r,r')F x (r')dr' = XF x (r) 



with the boundary conditions F\(0) = F X (R) = 0. If all the eigenvalues are negative, then 
5 2 S < and the critical point is a local entropy maximum. If one eigenvalue is positive, the 
critical point is an unstable saddle point. The point of marginal stability A c is determined by 
the condition that the largest eigenvalue is equal to zero (A = 0). We thus have to solve the 
differential equation 



(85) 
with 
(86) 

and F(0) 
rewritten 

(87) 



k d 



( 



d 



m dr \Aiipr 2 dr J Tr 



+ 



G 



F(r) 



V— (r) 



Mmc A T'(x) dr 



V= / ^( r ')F(r')dr', 
Jo dr 



F(R) = 0. Introducing the dimensionless variables defined in section |2]J, it can be 



d fe* d 



d£ 
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with 



X = 



and F(0) = F(a) = 0. As shown by Padmanabhan (1989), the solutions of the differential 
equation (|87| ) can be expressed in terms of the solutions of the Emden equation ( [52] ) as 



(89) 



F(0 = X 



1 



-(e 3 e-^-eV) + eV 



We can check that this function satisfies the boundary conditions F(0) = F(a) = 0. The point 
of marginal stability is obtained by substituting the solution ( j39|) in Eq. (|88|) . The integrations 
can be carried out and the solutions expressed in terms of the Milne variables u and v (see 
Padmanabhan (1989) for more details). In our semi-relativistic treatment, we obtain 



2u + u v - 7u + 3 - v Q F' '(fj,v )(u - 1) = 



(90) 

where we have used Eqs. (|64] ) and flT5"D to simplify the last term. 

In the classical limit fi — > +oo we recover the result of Padmanabhan (1989) 



(91) 

and in the formal limit \i 
(92) 



Au + 2uqVq _ ll^o + 3 
0, we find 

2u + v - 4 = 0. 



0. 



The intersections between the curve (^) and the spiral in the (u, v) plane determine the 
values of a for which a new mode of stability is lost (i.e., a new eigenvalue A becomes positive). 
Since the curve ( |90D passes through the singular sphere (u s ,v s ) = (1,2) at the center of the 
spiral, there is an infinity of intersections. The first intersection (for which a is minimum) 
corresponds to the point of marginal stability. We can show that the points determined by Eq. 
( |90| ) are precisely those for which A is extremum in agreement with the turning point analysis 
of Katz (1978). Indeed, differentiating the expression ([T(J) for A with respect to a, we get 



(93) 



dA 

da 



— = -[i T {/mvq)- 



da 



3 dvo 
Vq da 



1 duo uq dvo 
vq da Vq da 



Using equations ([72]) -( [73|) , we obtain 



(94) 



dA 

da 



avo 



2u + u v - 7u + 3 - fi v F' (iiv )(u - 1 



and we check that the condition ^ = is equivalent to Eq. (|90|) . 



da 



It is also easy to determine the form of the perturbation that triggers the instability at the 
critical point. According to Eq. ([ST|), the eigenfunction associated with the eigenvalue A = 
can be written 



(95) 



dp 
Po 



1 dF 
A^~dJ' 



where F(£) is given by Eq. (|89|) . Simplifying the derivative with the aid of Eq. 
express the perturbation profile in terms of the Milne variables ([TT]) as 



we can 



(96) 



Sp 
P 
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Figure 6: Graphical construction to determine the nodes of the perturbation profile 5p at the 
point of marginal stability. The construction is done explicitly for \i = 0, for which ct\ = 16.0 
and u = 0.817. There is only one zero satisfying £i < ct\ so that the perturbation profile does 
not present a "core-halo" structure. This property is maintained until \i > 1.61. 
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Figure 7: Density perturbation profile at the point of marginal stability as a function of the 
relativistic parameter. 
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The qualitative behavior of the perturbation profile can be studied without numerical inte- 
gration by a graphical construction (see Padmanabhan 1989). The density perturbation Sp 
becomes zero at the point (s) & such that uq = 3 — v(£i). On Fig. §, we first draw the line 
u = 3 — v . This lines passes through the singular sphere (u s = l,v s = 2) and also cuts the 
spiral at the points of vertical tangent (see Eq. (|72|)). In particular, the first intersection cor- 
responds to a* = 22.5 and (w*, u*) = (0.793,2.208). Then, we draw the line u = u = u(ai). 
The intersection between these two lines determines f(£i). The intersection between v = 
and the spiral determines the zeros of 5 p. For a > a*, there are two intersections satisfying 
£i < q;i so that the perturbation profile presents a "core-halo" structure. This is the case in 
particular in the classical limit p — > +oo for which «i = 34.4 (see Padmanabhan 1989). By 
contrast, for a < a*, there is only one intersection satisfying £i < a% so that the perturbation 
profile does not present a "core-halo" structure. This is the case in particular in the (formal) 
limit p = for which a\ = 16.0. The "core-halo" structure disappears for a± = a* = 22.5 
corresponding to a relativistic parameter p* = 1.61. Since our study is valid for p 1, we 
deduce that the density profile always presents a "core-halo" structure in the cases of physical 
interest. However, relativistic effects have the tendency to reduce the extent of the halo (see 
Fig. 0). 

The previous results are valid in the microcanonical ensemble in which the energy is fixed. 
In the canonical ensemble, we must consider maxima of the free energy J = S — (3E at fixed 
temperature. In that case, the condition of stability is given by Eq. fl8"5] ) with V = 0. Since the 
relativistic function F{x) does not appear anymore in the equations, we conclude that special 
relativity does not change the classical results in the canonical ensemble. In particular, the 
perturbation profile does not present a "core-halo" structure at the critical point (Chavanis 
2001). 



3 Isothermal gaseous spheres in general relativity 
3.1 The equations governing equilibrium 

We now address the structure and the stability of isothermal gas spheres in the context of 
general relativity. The Einstein field equations of general relativity are expressed as 

(97) R/J.U ~^9ftvR ^ Tftv, 

where R^ v is the Ricci tensor, T^ v the energy-momentum tensor and the metric tensor 
defined by 

(98) ds 2 = -g^dafdx", 

where ds is the invariant interval between two neighbouring space-time events. 

In the following, we shall restrict ourselves to spherically symmetric systems with motions, 
if any, only in the radial directions. Under these assumptions, the metric can be written in the 
form 

(99) ds 2 = e u dr 2 - r 2 {dd 2 + sin 2 Odtf) - e x dr 2 , t = ct, 

where v and A are functions of r and r only. The energy-momentum tensor is assumed to be 
that for a perfect fluid 

(100) = pg^ + {p + e)uV, 
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where = dx^/ds is the fluid four- velocity, p is the isotropic pressure and e is the energy 
density including the rest mass. 

The equations of general relativity governing the hydrostatic equilibrium of a spherical 
distribution of matter are well known. They are given by (see, e.g., Weinberg 1972) 

d , \ s 8irG 9 

101 re -A =1 _ r 2 e> 

dr c 4 

/,™n dp 1, .dv 

dr 2 dr 

e~ x dv 1 ,„ x. 8irG 
(103) __ = _ (1 _ e ->) + _ p . 



These equations can be deduced from the Einstein equations fl9~7|) . However, Eq. ( |102|) can 
be obtained more directly from the local law of energy-momentum conservation, D^T^ = 0, 
which is also contained in the Einstein equations. 
Equations (|101|) - fll03|) can be combined to give 

(,04) { 1 _^}^.4 (£ + P ){£^) + ^ pr 

[ c r J dr c [ r z c 

and 

(105) M(r) = — / er 2 dr. 



o 



These equations are known as the Oppenheimer-Volkoff (1939) equations. They extend the 
classical condition of hydrostatic equilibrium for a star in general relativity. Using equations 
(|101|) and ( |103| ), we find that the functions A(r) and z/(r) satisfy the relations 

2G 

(106) e~ x = 1 -M(r), 

rc 

and 

dv 1 + 4irpr 3 /M(r)c 2 



(107) 



dr r(rc 2 /2GM(r) - V 



In the empty space outside the star, p = e = 0. Therefore, if M = M(R) denotes the total 
mass-energy of the star, Eqs. (|106|) ( |1071 ) become for r > R 



_ A _ 2GM dv _ 1 
{ } 6 ~ rc 2 ' dr~ r(rc 2 /2GM - 1) ' 

The second equation is readily integrated in 

(109) z/ = ln(l-2GM/rc 2 ). 

Substituting the foregoing expressions for A and v in Eq. (0), we obtain the well-known 
Schwarzschild's form of the metric outside a star 

(no ^ = (i - - w + ^) - r^j-. 



19 



This metric is singular at 

, N 2GM 
HI r = ^^ = R s , 

cr 

where Rs is the Schwazschild radius appropriate to the mass M. This does not mean that 
spacetime is singular at that radius but only that this particular metric is. Indeed, the singu- 
larity can be removed by a judicious change of coordinate system (see, e.g., Weinberg 1972). 
When Rs > R, the star is a black hole and no particle or even light can leave the region 
R < r < Rs- However, in our case the discussion does not arise because Rs < R. Indeed, for a 
gaseous sphere in hydrostatic equilibrium, it can be shown that the radius of the configuration 
is necessarily restricted by the inequality (Buchdahl 1959) 

8 cr 8 

Therefore, the points exterior to the star always satisfy r > Rs- 
3.2 The equation of state 

To close the system of equations (|104|) - (|105|) , we need to specify an equation of state relating 
the pressure p to the energy density e. Quite generally, the first law of thermodynamics can be 
expressed as 

(113) d(^j = -Pd(^j + Tds, 
or, equivalently, 

(114) de = ^ dn + nTds, 

n 

where n is the baryon number density and s the entropy per baryon in rest frame. It must 
be completed by two equations of state p = p(n, s) and T = T(n, s). Then, Eq. ( |114j ) can be 
integrated to give e(n,s) (see, e.g., Misner et al. 1973). 



We shall assume in the following that the term Tds in Eq. flll4j ) can be neglected. This 
simplification arises in two different situations. In the case of neutron stars or white dwarfs, 
the thermal energy kT is much smaller than the Fermi energy, so the neutrons or the elec- 
trons are completely degenerate (kT <C Ef ermi ). On the other hand, in supermassive stars, 
temperature and entropy are important but convection keeps the star stired up and produces 
a uniform entropy distribution (s = Cst.). In these two important situations, the first law of 
thermodynamics reduces to 

(115) de = - 6 dn, 

n 

and we just require one equation of state p = p(n). We shall consider an equation of state of 
the form 

(116) p = Kn 1 , 



where K and 7 are constant. It is easy to check that the general solution of Eq. ( |1 1 5|) with 
Eq. (|11D is 

(117) e = Ap lh + : -p, 

7 - 1 



20 



where A is a constant. Systems obeying the pressure-energy density relation ( |117| ) with A ^ 
have been considered by Tooper (1965). Since the relation between e and p is essentially a 
power-law, these systems generalize the polytropes of Newtonian theory. We shall here consider 
the limiting situation A = 0. In that case, Eq. ( |117j ) reduces to the so-called "gamma law" 
equation of state 

(118) p = qe with q = 7 — 1. 

Since the relation between p and e is linear, this equation of state extends the theory of isother- 
mal spheres to the context of general relativity. As noted by Chandrasekhar (1972), for this 
equation of state, the Oppenheimer-Volkoff equations become mathematically similar to those 
describing a classical isothermal gas (the so-called Emden equations). 

An equation of state of the form (|118|) has been introduced in different situations: 
(i) This equation of state prevails in the highly energetic core of neutron stars where the 
matter is ultra-relativistic and completely degenerate (Oppenheimer & Volkoff 1939, Misner 
& Zapolsky 1964, Meltzer & Thorne 1966). If the system is modelled as a pure collection of 
noninteracting fermions, standard theory leads to (Chandraskhar 1942) 

(119) p=U (7 = 4/3). 



In that case, the constant K which appears in Eq. (116) is explicitly given by K = |(-) ly/3 /ic, 



where h is the Planck constant. Other versions of the equation of state have attempted to take 
into account nucleon-nucleon interactions or a spectrum of baryon species. Equations of the 
form ( |118| ) but with numerical coefficient other than 1/3 have sometimes been suggested. For 
example Ambatsumyan & Saakyan (1961) have used an equation with q = 1/13 (7 = 14/13) 
and Tsuruta Sz Cameron (1966) with q — 1. The equation of state 

(120) p = e, 

was also introduced by Zel'dovich (1962) for a gas of baryons interacting through a vector 
meson field. This represents an extremely high pressure since the sound velocity, (dp/deY^ 2 c, 
is equal to the velocity of light for this value of q. This is clearly an upper bound. In the 



following, we shall consider a general equation of state of the form ( 118 ) with q in the range 
< q < 1. 

(ii) Bisnovatyi-Kogan & Zel'dovich (1969) and Bisnovatyi-Kogan & Thorne (1970) have 
considered the general relativistic equilibrium configuration of a gas whose temperature T is 
constant everywhere. In that case, the equation of state is given by (see section □) 



'121) P=— T 

Tt K 3(t^) - A 2 (^ 



which is again of the form ( 118 ). In particular, q = \ when T — > +00. However, as Bisnovatyi 



Kogan & Thorne emphasize, the equation of state ( |121[ ) can only describe a "local" thermo- 
dynamical equilibrium. Indeed, "global" thermodynamical equilibrium in general relativity 
requires that the redshifted temperature e v ^' 2 T be uniform throughout the medium (Tolman 
1934). Thus, the spheres described by Eq. (|121 ) are not "isothermal" in the general- relativistic 



sense; heat will slowly diffuse inward in them, upsetting the condition T = Cst. and trying 
to establish e v( ^^ 2 T{r) = Cst. instead. However, since the hydrostatic equations ( |104[ ) ( |105l ) 
with the equation of state (|118 ) are mathematically similar to the classical Emden equation 
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for isothermal gas spheres, we shall call the equation of state (|118| ) "isothermal" , following the 
terminology of Chandrasekhar (1972), although this is only correct in a local sense. 

(iii) The equation of state (|118[) was also proposed by Saslaw et al. (1996) in a cosmological 
context to model the ultimate state of an Einstein-de Sitter universe that undergoes a phase 
transition caused by gravitational clustering. This phase transition can lead to the growth of a 
centrally concentrated distribution of matter so that the universe would pass from a statistically 
homogeneous state to a state of rotational symmetry around one point only. The resulting 
configuration neither expands nor contracts, so the global solution is stationary. According 
to Saslaw et al. (1996), this "isothermal universe" would represent the ultimate astrophysical 
prediction. 

(iv) Note finally that the equation of state p = |e prevails in a medium where the pressure 
is entirely due to radiation (Chandrasekhar 1942). 

3.3 The general relativist ic Emden equation 

Considering the equation of state ( |118| ), we shall introduce the dimensionless variables £, if) and 
M(£) by the relations 

< i22) ^{«^W /2 «, 



and 



In terms of the variables ip and £, Eqs. ( |104|) and ( |105| ) can be reduced to the following 
dimensionless forms (Chandrasekhar 1972) 



and 



(125) ^-^ = ee 



In addition, the metric functions determined by equations (|106|) and (|102|) can be expressed as 



(126) e-^-^fl, ^JS_ m + C s,, 

1 + q 4 1 + 9 

where the constant is determined by the matching with the outer Schwarzschild solution 
at r = R. 

The Newtonian limit corresponds to q —>■ 0. In that limit, Eqs. (|124j) -( 125 ) reduce to 



and they combine to give Emden's equation 
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Therefore, Eqs.( |124|) - (|125|) represent the general relativistic equivalent of the Emden equation. 
Like in the Newtonian case, they admit a singular solution of the form (Chandrasekhar 1972) 

(129) <^' = % where Q 2(1 + q) 



The metric associated with the singular isothermal sphere is given explicitly by 



(130) e u = A£—* } e A = l 



where A is an unimportant constant. Considering now the regular solutions of Eqs. ( |124| )-( |125|) , 
we can always suppose that e represents the energy density at the center of the configuration. 
Then, Eqs. (|124j) - (|125|) must be solved with the boundary conditions 

(131) ^(0) = V'(0) = 0. 

These solutions must be computed numerically. However, it is possible to determine their 
asymptotic behaviors explicitly. For £ — > 0, 

(132) ij = a£ 2 + b^ + ..., 
with 

(133) „ = !±* t= _W-2» + 3)(l + 3,) 
V ; 6 ' 360(1 +q) 

and for £ — > +oo (Chandrasekhar 1972) 



(134) ^=^l + ^ w(1+g) °°s( ; 2(1 ' +g )' + * 

The curve ( |134|) intersects the singular solution ( |129| ) infinitely often at points that asymp- 
totically increase geometrically in the ratio 1 : ex P (7+42^2)1/2 ( see Fig- ||)- Since the energy 
density e falls off as £ -2 at large distances, the total mass-energy M is infinite. In practice, 
this "infinite mass problem" does not arise because the isothermal equation of state (|118|) only 
holds in a finite region of space (e.g., the inner regions of a neutron star). For simplicity, we 
shall remedy this difficulty by assuming that our system is confined within a box of radius 
R, like in section |2|. This is clearly an idealization but it provides a well-posed model which 
captures the essential features of the system (see below) and which can be be studied in great 
detail without any further approximation. 

If the system is enclosed within a box, the solutions of Eqs. (|124|) ( |125| ) must be terminated 
at different radii a given by 

(135) a= {j^±j)} 1/2 fl . 

It should be noted that a is a measure of the central energy density e . Instead of a, we might 
prefer to consider the density contrast 

(136) n = -^ = e^ a) . 
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Figure 8: Density profiles of the general relativistic Emden equation for different values of q. 
The profiles behaves like £ -2 at large distances. 



The density contrast is a monotonic function of a varying from 1Z — 1 (homogeneous system, 
a = 0) to 1Z — > +oo (singular sphere, a — > +oo). It is sometimes of interest to express the 
results in terms of the redshift (see, e.g., Weinberg 1972) 

(137) 2 =y = Wr)rl/2 - 1 = e " Kr>/2 - 1 

Using Eq. ( |12(j| ) with the boundary condition ( |109| ) at r = R, the redshift of a spectral line 
emitted from an isothermal sphere is given by 

(138) ^expl^-M^-m)]-!- 



l + q 



3.4 The Milne variables 



As in the Newtonian theory, it will be convenient in the following to introduce the Milne 
variables 

(139) U =%T ' and v = ^'- 



In terms of these variables, the system of equations (|124|)(|125|) can be reduced to a single first 
order differential equation (Chandrasekhar 1972) 



udv -l-^»+(l + 3 t )«+'^ m) + ^ 
vdu 3 - - (1 + 3q)u - ^uv - ^uv* 

For ^ — > 0, one has 

(141) u= JL-(£ + !)«. + ..., „ = 2 ^ + ..., 



(140) 
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and for £ — > +00 
(142) 



u — > ti, = — 



2 ' 



V — > u, 



3 




Figure 9: The (w, t>) plane for isothermal gaseous spheres in general relativity and for different 
values of q. The value q = corresponds to the Newtonian limit. 

The solution curve in the (u, v) plane is parametrized by £. Starting from the point (u, v) = 
1 " 0) for £ = with a slope (dv/du) = — |^r, the solution curve spirals indefinitely around 



the point (u s , u a ) = (§,2), corresponding to the singular sphere, as £ — > +00. All isothermal 
spheres must necessary lie on this curve (see Fig. ||). For bounded isothermal spheres, £ must 
be terminated at the box radius a. Clearly, the spiral behavior of the (u, v) curve can be 
ascribed to the oscillating behavior of the solution (|134j) as £ — > +00. An explicit equation for 
the spiral (valid for £ — > +00) can be obtained by substituting the asymptotic expansion ( |134j) 
in the Milne variables ( p.39|) like in the Newtonian case (see Chavanis 2001). 

3.5 Oscillatory behavior of the Mass-density profile 

According to Eq. (|123[) , the relation between the total mass M of the configuration and the 
central energy density eo (through the parameter a) is given by 

/ \ q M(a)Rc 2 

143 M=— 

y 1 1+q a G 

Solving for M(£) is Eq. (|T2§, we get 

(M4) M(a) _ afja) - ga 2 e-^ ( 2q 



a l+paip'(a) \ l + q, 

This relation can be expressed very simply in terms of the values of the Milne variables u and 
v at the normalized box radius a. Writing u = u(a) and v (a) and using Eq. ( |139|) , we obtain 

2GM _ pv (l-qu ) 
1 } X ~ Rc 2 l+pv ' 
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Figure 10: Total mass-energy x versus central energy density a for q = 1/3. 



The curve x( a ) starts from x — for a = and oscillates around its asymptotic value 
Xs = pQ = n+nV+dn (corresponding to the singular sphere) as a — > +oo (see Fig. [10]). This 



(l+<?) 2 +4g 

relation between the total mass and the central density is very similar to the corresponding 
one for neutron star models (Oppenheimer & Volkoff 1939, Misner & Zapolsky 1964) and dense 
stellar clusters at statistical equilibrium (Bisnovatyi et al. 1998). In particular, there exists 
mass peaks in the diagram. In fact, this oscillatory behavior should not cause surprise since 
the same phenomenon exists for a classical isothermal gas [see in particular Fig. 7 of Chavanis 
(2001) giving the mass-density relation for a fixed temperature]. 
From Eqs. (PgpflH^) and Ql25|) , one has 



(146) 



dx _ 
da 



M'(a) 



v- 



M(a) 



pae 



a a* 

In terms of the Milne variables, it can be rewritten 



V 



1 + paijj'(a) 



(147) 



dx 
da 



u v 



v (1- qu ) 



l+pv 

Therefore, the extrema of the curve x{ a )i determined by the condition dx/da 



0, satisfy 



(148) 



pv 



--q-1. 

u 



This equation defines a hyperbole in the 



. u. v 



plane. For u — * + , v ~ — > +oo and for 



+ 00, Vq 



(<?+l) 2 
2q 



The intersections between this curve and the spiral (see Fig. [TT|) 



determines the values of a for which x is an extremum. Since the curve ( |148| ) passes through 
the singular sphere (u s ,v s ), there is an infinity of intersections with the spiral, resulting in an 
infinity of oscillations in Fig. [T0| . 

The first value ol\ corresponds to the global maximum of x- Therefore, isothermal spheres 
in general relativity exist only provided that 

pu(ai)(l - gtt(ai)) 



(149) 



X < Xc 



1 +pv(a 1 ) 
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Figure 11: Graphical construction to determine the critical density a± for q = 1/3. 

This implies in particular that, for a given radius, confined isothermal gas spheres exist only 
below a limiting mass 

Rc 2 

(150) M < M, - \ 



2G 



This result is of course related to the existence of a limiting mass for neutron stars as evidenced 
by Oppenheimer & Volkoff (1939). In their study, the value of the radius R is determined by a 
proper modelling of the envelope. However, our results indicate that the essential properties of 
neutron stars: oscillatory behavior of the mass-density profile, limiting mass and spiral behavior 
of the mass-radius diagram (see section |3.6| ) are due primarily to their isothermal core and not 
to their envelope. Similar observations have been made by Yabushita (1974) who considered 
the case of an isothermal gaseous sphere surrounded by a medium exerting on it a constant 
pressure. 

The critical parameter Xc, the critical density contrast 1Z C and the critical redshifts z c 
emitted from the center and from the boundary of the isothermal configuration are ploted as a 
function of q on Figs. |T2[|T^4TJ. The corresponding values of a\ are indicated on Fig. |TB] (full 



line). For the Newtonian case (q = 0), we recover the classical values ocx = 8.99 and 7Z C = 32.1 
obtained in the canonical ensemble. When q is increased, the critical density a± is lowered 
so that instability occurs sooner than in the Newtonian case. It should be noted, however, 
that the critical density contrast 1Z C is not a monotonous function of q. This is due to the 
deformation of the spiral in the (u, v) plane when we vary the relativistic parameter q. We find 
that the value of q for which 1Z C is minimum is close to 1/3, the typical value corresponding 
to neutron stars. However, the variation of the critical density contrast with q is not very 
important (1Z C ~ 20 — 30 in the whole range of parameters) so that the critical central redshift 

(151) z° = (1 - XcY XI2 TlP ~ 1, 

obtained from Eqs. ( |138| ) ( |136|) ( |145|) , increases monotonically with q. 



27 




Figure 12: Critical parameter Xc as a function of q. The classical limit q — > is discussed in 
section [O 
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Figure 14: Critical redshifts z c emitted from the center and from the boundary of an isothermal 
configuration as a function of q. 

3.6 The mass-radius diagram 

It is well-known that, for high central densities, the mass-radius diagram for neutron stars 
presents a spiral behavior. We show in this section that we can reproduce this behavior within 
our simple "box" model. To that purpose, we shall consider configurations with different masses 
and radii but with the same density e(R) at the boundary of the domain. According to Eqs. 
(|135|) and ( |136| ), we can express the radius R as a function of the parameter a by the relation 



(152) 4- = «e"^ (Q)/2 , 

Ro 

where we have introduced a typical radius 

A 1/2 



(153) R 



qc 



47rGe(R)(l + q) 



On the other hand, using Eq. (|145|) , the dimensionless mass appropriate to the present context 
is given by 

(154) ^ = X (a)^(a), 

1V1 Ho 

where 

(155) M '" Ri) 



2G 



Eqs. ( |152j ) ( |1 54] ) determine the mass-radius relation for an isothermal gas with a fixed energy 
density at its boundary. The M—R curve, represented on Fig. [H| is parametrized by the central 
density a. For comparison with real objects (e.g., neutron stars) we note that the assumption 
of an isothermal core is valid only for sufficiently large central densities, typically from the 
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Figure 15: The mass-radius diagram for an isothermal gas with a fixed energy density at the 
boundary (e(R) = 10 15 gc 2 /cm 3 ). We have indicated the value of the central density eo at the 
turning points. 



point encircled. From that point, our model reproduces qualitatively the spiral behavior of the 
mass-radius diagram for neutron stars. Using the asymptotic expansion ( |134j ) of the function 
t/>(£)> ^ is possible to determine an explicit equation for the spiral (valid for high densities) in 
the form 

( 156 ) ir Ql/2 { 1+ ^^ cos 



V 2(1 + q) 



M „,,,f. A 



P Q 3/2 {l + 



M ™ [ 4(l + 5g)a( 1+3 9)/ 2 ( 1 +9) 



[l + q)(7 + 42q-q 2 ) 



2\l/2 



. f(7 + A2q-q 2 ) 1 /\ \ . in 2n / (7 + A2q - q 2 ) 1 ' 2 , 

(157) 

In particular, the terminal point of the spiral, corresponding to the singular solution (|129 ), is 
given by 

(158) £ = Q 1/2 , § = pQ"\ 

The previous results are valid for any value of the parameter q. We shall now specialize on the 
case of neutron stars for which q = 1/3. For estimating e(R), we shall adopt the typical value 
given by Misner & Zapolsky (1964) 

(159) e(R) = 10 15 gc 2 /cm 3 , 

which corresponds in an ideal neutron gas to a Fermi energy of 130 MeV. This is the typical 
value of the energy density in the regime where the kinetic energy of the neutrons is of the same 
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order as their rest mass. This regime corresponds to the transition between the isothermal core 
and the envelope and it is therefore relevant to adopt this value at r = R. For this density, 



(160) Rq = 5.2 km, M = 1.8 M & . 

From the diagram of Fig. [15|, we find that the maximum mass and the corresponding radius 



have the values 

(161) R c = l.LRo, M c = 0.5M . 

The value of the maximum mass in our model M c ~ O.9M is close to the more exact value 
~ O.7M found by Oppenheimer & Volkoff (1939) by introducing a more general equation of 
state valid also in the envelope. The value of the corresponding radius R c = 6 km is slightly 
smaller than the value ~ 9.6 km usually reported (see, e.g., Weinberg 1972) and the difference 
in attributed in part to the presence of the envelope which surrounds the isothermal region of 
size R c . According to Eq. (|138|) , the fractional redshift of a spectral line emitted from the edge 



of the isothermal region (for the configuration of maximum mass) is 
(162) ^ )= ( 1 -^ £ )" /2 - 1 = - 35 - 

and the corresponding central redshift Zq = 1.4. Finally, we find that the central density at 
the mass peak is = 10 16 gc 2 /cm 3 is good agreement with the reference model of Meltzer 
& Thorne (1966) giving a comparable value. Other central densitites are indicated on Fig. 



15 . Our curve matches relatively well (qualitatively and semi-quantitatively) the diagram of 



Meltzer & Thorne (1966) for central densities eo > 4 10 15 gc 2 /cm 3 . 



4 Dynamical instability of isothermal gas spheres in gen- 
eral relativity 

4.1 The condition of marginal stability 

The stability of relativistic stars against radial perturbations is usually studied with the equa- 
tion of pulsations derived by Chandrasekhar (1964). Then, the stability of stars can be reduced 
to a Sturm-Liouville problem, which is usually solved numerically (see, e.g., Bardeen et al. 
1966). The method is quite general but it is not the most efficient in the present situation. 
Indeed, by formulating the problem differently, using the form of the equation of pulsations 
given by Yabushita (1973) and making use of the Milne variables, it is possible to study the 
stability of bounded isothermal spheres analytically, like in the Newtonian case. 

Let us consider a small perturbation 5e around a configuration of isothermal equilibrium 
and let us introduce the mass perturbation /(£) within the sphere of radius £ such that 

163 _ = T~C2lc' 

eo 4tt£ 2 d£ 

Yabushita (1973,1974) has shown that the equation for radial pulsations satisfied by the function 
/ can be written 



i + 



2g£ #1 £ a 2 c* x _ v 



f = e x - u f 

J 4nGe (l + q) J ' 
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where a is the period of oscillations defined by 5e ~ e act . Clearly, the condition of instability 
corresponds to o 2 > 0. On the other hand, the velocity profile of the perturbation is given by 

(165) Su 1 = - — -LL e -»/ 2 . 

q + 1 47re r l 

If the system is confined within a box, we must require that Su 1 = at r = R. The function 
/(£) must therefore satisfy the boundary conditions 

(166) /(0) = f(a) = 0. 

These boundary conditions are equivalent to the conservation of mass-energy M. In the New- 
tonian limit q —>■ 0, the equation of pulsations ( 164 ) reduces to the one considered by Yabushita 



(1968) and Chavanis (2001). We shall be particularly interested by the point of marginal sta- 
bility (a = 0) in the series of equilibrium at which the system becomes unstable. We then have 
to solve the equation 



l + qdi 



0. 



with F(0) = F(a) = 0. If we denote by £ the differential operator that appears in Eq. (|167D , 
it is possible to show that (Yabushita 1974) 

Therefore, the eigenfunction associated with a = is 
(169) F(0 = Cl (e 3 e^-M(0). 

In the classical limit, we recover the result £ 3 e - ^ — £ 2 ip' derived from the analysis of the 
second order variations of the free energy (in the canonical ensemble) or from the Navier-Stokes 
equations (Chavanis 2001). The condition -F(O) = is automatically satisfied. The condition 
F(a) = determines the values of a at which a new mode of stability is lost. According to Eq. 
flU), we get 



(170) M(a) = a s e 



3 



Using Eq. ( |144| ) and introducing the Milne variables ( |139| ), it is straightforward to check that 
this relation is equivalent to 

(171) pv Q = q-1, 

which is precisely Eq. ( |148| ). Therefore, a new eigenvalue a 2 becomes positive at the points 
where \ is extremum. In particular, the first instability sets in precisely at the point of maxi- 
mum mass. Isothermal configurations with a central density a > a±, a density contrast TZ > 1Z C 
or a central redshift z > Zq are unstable with respect to radial perturbations. Similar results 
were obtained by Misner & Zapolsky (1964) who numerically solved the Chandrasekhar equa- 
tion for radial pulsations applied to their neutron star model. Their method was, however, 
limited to the fundamental mode and they could not demonstrate that secondary mass peaks 
correspond to the change of sign for the next eigenvalues (see their discussion). Our method 
gives an unambiguous definite answer to that question, at least in the framework of our model. 
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It has to be noted that these results differ from those obtained by Yabushita (1974) who found 
that instability occurs before the first mass peak (except in the Newtonian case q = 0). Now, 
the three models only differ in the boundary conditions: (i) Yabushita's gas spheres have en- 
velopes which exert a constant pressure, (ii) In our model, the volume of the isothermal region 
is assumed to be fixed, (iii) Misner & Zapolsky introduce a more realistic envelope in which 
the equation of state has a physically admissible form. Our results suggest that imposing a 
fixed volume (Antonov point of view) is a better approximation than imposing a fixed pressure 
(Bonnor point of view) since it provides the same condition of stability as the reference model 
of Misner & Zapolsky (1964) and other works. 



4.2 Perturbation profiles and secondary instabilities 



The values of a at which a new mode of instability occurs are determined by Eq. ( |171|) . For 
large values of a (high order modes), we can introduce the asymptotic expansion ( |134j ) in 
the Milne variables. Then, after some algebra, we find that the values of a n at which a new 
instability occurs are asymptotically given by 



(7 + 42^-gY/2 / g 2 + 2Qg + 3 
(172) — — - in a n + o = arctan ( — ^7^— — — ON1 /n | + mr. 



- q 2 ) 1 ' 2 

2{lTq) ' " ^™V(l + g)(7 + 42g-g 2 )V2 

In continuity with the Newtonian limit, they follow the geometric progression 

27r(l + 9 ) 



(173) 



g(7+42g- ? 2) 1 /2 



in 



-00, integer). 



For q = 0, 1/3, 1 the ratio in the geometric progression is respectively 10.74, 6.25 and 6.13. It 
is minimum for q = 7/11 with a value 5.93. 

The profile of the energy density perturbation associated with each instability is obtained 
by substituting Eq. ( |169j ) in Eq. ( |163| ). Simplifying the derivative with the aid of Eq. ( |125| ) 
and introducing the Milne variable v, we obtain 



(174) 



5e 



47T 



'2-v) 



The zeros of the profile 5e are determined by the intersections between the spiral in the (u, v) 
plane and the curve v — 2 (see Fig. |ll|). The first two zeros are represented on Fig. [16] 
for different values of q. We see that £2 is always larger than a\ so it does not satisfy the 
requirement < £ < a. Therefore, the first mode of instability has only one node at £ = £1. 
In other words, the perturbation profile that triggers the instability of an isothermal sphere in 
general relativity does not present a "core-halo" structure (see Fig. 0) in continuity with the 
Newtonian study in the canonical ensemble (Chavanis 2001). 

The high order modes of instability have several nodes. Their positions (for large £j) are 
obtained by expanding the Milne variable v with the aid of Eq. (|134 ) and substituting the 
resulting expression in the condition v — 2. This yields 



(175) 



(7 + A2q-q 



In + 5 = — arctan 



- 



1 + 3g 



+ iir. 



2(1 + q) ' " — V(7 + 42g-g 2 )V2 

Therefore, the zeros of the perturbation profiles also follow a geometric progression 



(176) 



6 



2?r(l + g) 



g (7+42q-(^ 



+00, integer), 
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Figure 16: Critical parameter a.\ as a function of q (full line) together with the first two zeros 
of the perturbation profile (|174j) satisfying v(£i) = 2 (dash lines). Since & > &i, this solution 
must be rejected. 
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Figure 17: Perturbation profiles at the point of marginal stability for different values of q. The 
perturbation profile does not present a core-halo structure. 



34 



with the same ratio as the dimensionless radii a n . 

We can also determine the velocity profile (165) associated with each mode of instabil- 



ity. Introducing the velocity of sound c s = q x l 2 c and expressing the result in terms of the 
dimensionless variables defined in section pO, we get 



. , Su 1 co' tm. „ tJ . s 

(177) --4^7W el+ ' m 

where F(£) is given by Eq. ( |169|) at the critical points. In writing Eq. (|177|) , it is implicitely 
understood that we are just at the onset of instability {&' = + ), so that Eq. ( |177| ) is applicable 
with a' > (the velocity profile at the critical points, i.e. a' = 0, is simply 5u l = 0). In terms 
of the Milne variables, the velocity profile ( |177j) can be written 



,,„„, Su 1 ca' ,, Ml u ( 1 

178 = - , , —r T Cii/e^j—- [pv-- + q + l 

c s 47r(g + 1) 1 +pv \ u 

For q = 0, we recover the expression obtained in the Newtonian context from the Navier-Stokes 
equations (Chavanis 2001). Since F(0) = F(a) = 0, the velocity always vanishes at the center 
and at the boundary of the domain. The other zeros, for which Ffe) = 0, correspond precisely 
to the values of a at which a new mode of instability occurs, i.e. ^ = on- The n th mode of 
instability therefore has n+1 zeros £ = 0, £ = cti, £ = ct 2 ,..., £ = a n . For high order modes of 
instability, the zeros follow asymptotically the geometric progression ([173|). 



4.3 The baryon number: another interpretation of the instability 

It has often been suggested on the basis of phenomenological arguments that the baryon number 
N, or equivalently the rest mass Mq = Nm, is also maximum at the onset of instability. We 
wish to check this prediction explicitly in the case of an isothermal gas sphere. 

The baryon number is obtained by multiplying the baryon number density n by the proper 
volume element e A( - r;//2 ] 4irr 2 'dr and integrating the resulting expression over the whole configu- 
ration. Using Eq. (|106[) , we get (see, e.g., Weinberg 1972) 



(179) 



N 



n 



2GM(r) 



rc 



-1/2 



4:7ir 2 dr. 



Now, according to the equations of state (|116D - (|118|) , the baryon number density is related to 
the energy density by 



180) 



n ' 



K 



Introducing the dimensionless variables defined in section |3.3|, we find that the baryon number 
can be written 



181) 



N 



3q+l 

a l +i Jo 



e !+« 



with the normalization constant 
(182) N, 
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Figure 18: Baryon number and total mass-energy as a function of the central density (through 
the parameter a) for q = 1/3. The peaks occur for the same values of the central density. 
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The baryon number (|181| ) is ploted versus the central density on Fig. [18] for q = 1/3. 



It exhibits damped oscillations, like the total mass energy, and tends asymptotically to the 
value -^iQ^i^ — pQ) 1 ^ 2 as a — ► +00 (singular sphere). It is observed numerically that the 
peaks of mass-energy and of baryon number occur for the same values of the central density 
(we have checked this property for various values of q). This is also clear from Fig. [19] where 
the mass-energy M is ploted as a function of the baryon number N. From this figure, we can 
understand physically why the system becomes unstable after the first mass peak. Indeed, for a 
given rest mass Nmc 2 there exist several solutions with different energy Mc 2 . Only the solution 
with the lowest energy ( "fundamental" state) is stable with respect to radial perturbations (see 
Weinberg 1972). We expect therefore some phase transitions to occur for a state of high energy 



to a state of lower energy (see the arrow on Fig. 19). Therefore, only the lower branch is stable 
and it corresponds to the configurations prior to the first mass peak (we come to the same 
conclusion by minimizing the binding energy E = (M — M )c 2 at fixed mass M). These results 
are consistent with those found by Tooper (1964) for polytropic gas spheres. 

4.4 Characteristic length scales 

We shall now reformulate the previous results in a different form by introducing appropriate 
length scales. This will provide a direct comparison with our discussion in the Newtonian 
context (Chavanis 2001). In general relativity, a gaseous sphere of mass M is unstable if its 
radius is smaller than a multiple of the Schwarzschild radius (Chandrasekhar 1964). For a 
confined isothermal gas, this condition is explicitly given by Eq. ( |150| ) which can be rewritten 

, s 1 2GM 

183 R < —r^Rs, Rs = —3-. 

On the other hand, in Newtonian gravity, the condition of gravitational collapse is that the size 
of the domain be larger than the Jeans length Lj. For example, in the case of an isothermal gas 
confined within a box, the condition of instability reads (de Vega & Sanchez 2001, Chavanis 
2001) 

It should be emphasized that these inequalities characterize the absence of hydrostatic equi- 
librium for an isothermal sphere. This is more stringent than the instability of an equilibrium 
configuration (see below). Although the sense of these inequalities is opposite, there is no con- 
tradiction since the size of the domain enters explicitly in the definition of the Jeans length 
(through the average density). For q — » 0, the two criteria become equivalent as we now show. 
In the classical limit, the rest mass dominates the kinetic energy so that e = pc 2 where p = nm 
is the classical mass density. Therefore, the equation of state (|118|) becomes p = qpc 2 and it 



must be compared with the equation of state of a classical isothermal gas p = ^kT. This yields 

kT 1 

(185) q = — = — — ((3mc 2 - +00). 

mcr pmc z 



Now, in the classical limit c — > +00, the critical parameter 1/xc given by Eq. (|149|) , diverges 
like 1/xc — l/2gu(ai) — » +00 while the Schwarzschild radius Rs goes to zero in such a way 
that the product (l/xc)Rs remains finite. Therefore, the inequality ( |183j ) is reversed 



(3GMm 3R 3 v(a x 



(186) R < = ^R> 

v(ai) v(ai)Lj 
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and we recover the result ( |184j) with v(oti) = 2.52. It can be noted that the same criterion can 
be written 



(187) M >M C = 2.52^^. 

Gm 

Therefore, a classical gaseous sphere of radius R and temperature T can exist only below a 
limiting mass M c given by formula ( |187| ). Obviously, the physical origin of this limiting mass is 
the same as for the limiting mass of neutron stars but with a completely different interpretation 
of the parameters. 

It is also relevant to introduce another length scale 



(188) U = ( — 



\AirGe (l + q) 



1/2 



which gives a good estimate of the "core radius" of a relativistic isothermal gas sphere. In the 
classical limit, L c reduces to the King's length (Binney & Tremaine 1987). In terms of the core 
radius, the parameter a can be written 

(189) a = 3^- 

The condition of instability a > a% is therefore equivalent to 

(190) R > yL c . 

Above this threshold, it is still possible to construct isothermal gas spheres in hydrostatic 
equilibrium but they are unstable. Secondary instabilities occur for radii that follow a geometric 
progression 



(191) R n 



2n(l + q) 
e (7+42q- 9 2 ) 1 / 2 



L c (n — > +oo, integer). 



In the Newtonian context, this hierarchy of scales was interpreted by Semelin et al. (1999) and 
Chavanis (2001) in relation with the fragmentation of the interstellar medium and the possible 
origin of its fractal structure. 

The instability criteria ( |183|) and (|190|) , expressed in terms of distinct length scales, can also 



be expressed in terms of the redshift, with a similar distinction between the central redshift 
and the surface redshift. According to Eq. ( |162| ), the surface redshift z(R) is an increasing 
function of the mass energy \ while the central redshift z , deduced from Eq. ( |138| ), is an 
increasing function of the central density a. Therefore, the surface redshift must necessarily 
satisfy z(R) < z(R) c since no hydrostatic equilibrium of isothermal spheres can exist above this 
threshold while the central redshift must satisfy zq < Zq because configurations with zq > Zq 
exist mathematically but are unstable. 



5 Conclusion 

In this paper, we have studied the structure and the stability of isothermal gas spheres in the 
framework of general relativity. We have found that the relativistic isothermal spheres (like 
neutron cores) exhibit the same kind of behaviors as Newtonian isothermal spheres but with a 
different interpretation of the parameters. In this analogy, the critical energy and the critical 
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temperature found by Antonov (1962) and Lynden-Bell & Wood (1968) for an isothermal gas 
are the classical equivalent of the limiting mass for neutron stars discovered by Oppenheimer 
& Volkoff (1939) [there also exists a limiting mass for classical isothermal gas spheres, see 
Eq. ( |187| )]. Similarly, the spiral behavior of the temperature-energy curve for a classical gas 
has the same origin as the mass-radius diagram in neutron stars. This spiral behavior, as 
well as the damped oscillations of the mass-density profile, are general features of isothermal 
configurations. 

Although our model can reproduce qualitatively all the features of neutron stars, it is 
however limited in its applications since the isothermal core is surrounded by an artificial 
"box" instead of a more physical envelope like in the studies of Oppenheimer & Volkoff (1939), 
Misner & Zapolsky (1964) and Meltzer & Thorne (1966). However, this simplification allows 
us to study the stability problem analytically (or with graphical constructions) without any 
further approximation. This is a useful complement to the more elaborate works of the previous 
authors who had to solve the pulsation equation numerically or with approximations. This is 
also complementary to the study of Yabushita (1974) who considered isothermal gas spheres 
surrounded by an envelope exerting a constant pressure. The boundary conditions are of 
importance since Yabushita came to the conclusion that the onset of instability occurs before 
the first mass peak. In contrast, our model gives a stability criterion consistent with that of 
Misner & Zapolsky (1964) and other works. Furthermore, if we assume that the energy density 
at the edge of the isothermal core has a value ~ 10 15 gc 2 /cm 3 (a typical prediction of nuclear 
models), we can reproduce almost quantitatively the results usually reported for neutron stars. 

It should be stressed that the analogy between neutron stars and isothermal spheres is a pure 
effect of general relativity. When gravity is treated in the Newtonian framework, the classical 
condition of hydrostatic equilibrium requires a relationship between the pressure and the mass 
density p = ran which, for dense matter (degenerate and relativistic), is a power law with 
index 7 = 4/3 (see Chandrasekhar 1942). However, when gravity is treated in the framework 
of general relativity, the Oppenheimer- Volkoff equations require a relationship between the 
pressure and the mass-energy density e which, for dense matter, is linear. Therefore, the core of 
neutron stars (treated with general relativity) is "isothermal" while the core of white dwarf stars 
approaching the limiting mass (treated in a Newtonian framework) is "polytropic" , although 
the same equation of state is used (that for a completely degenerate and ultra-relativistic ideal 
Fermi gas). This is the intrinsic reason why the Mass- Radius relation for neutron stars exhibits 
a spiral behavior (Meltzner & Thorne 1966) while the Mass-Radius relation for white dwarf 
stars is monotonous (Chandrasekhar 1942). These remarks may increase the interest of studying 
isothermal gas spheres both in Newtonian mechanics and in general relativity. 
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